---
title: "conjoint analysis"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(tidyverse)
require(cjoint)
library(cregg)
library(scales)
library(patchwork)
library(extrafont)
require(memisc)

dat <- read_csv("data/cj_data.csv")
dat <- na.exclude(dat)

str(dat)
dat[,c(43:45)] <- data.frame(lapply(dat[,c(43:45)],as.factor))
dat[,c(47:48)] <- data.frame(lapply(dat[,c(47:48)],as.factor))

# Islamism
dat <- dat |> 
  mutate(Islamism = case_when(Q2_1 >= 67 ~ "High",
                              Q2_1 <= 66 & Q2_1 >= 34 ~ "Middle",
                              Q2_1 <= 33 ~ "Low"))
dat$Islamism <- as.factor(dat$Islamism)

# Liberalism
dat <- dat |> 
  mutate(Liberalism = case_when(Q2_2 >= 67 ~ "High",
                                Q2_2 <= 66 & Q2_2 >= 34 ~ "Middle",
                                Q2_2 <= 33 ~ "Low"))
dat$Liberalism <- as.factor(dat$Liberalism)

# Tribalism
dat <- dat |> 
  mutate(Tribalism = case_when(Q2_3 >= 67 ~ "High",
                               Q2_3 <= 66 & Q2_3 >= 34 ~ "Middle",
                               Q2_3 <= 33 ~ "Low"))
dat$Tribalism <- as.factor(dat$Tribalism)

# Child marriage
dat <- dat |> 
  mutate(Childe_marriage = case_when(Q14 >= 3 ~ "Agree",
                                     Q14 == 2 ~ "Neither",
                                     Q14 == 1 ~ "Disagree"))
dat$Childe_marriage <- as.factor(dat$Childe_marriage)

# Tribal/Religious
dat <- dat |> 
  mutate(Tribal_religious = case_when(Q15 == 1 ~ "Tribal",
                                      Q15 == 2 ~ "Religious"))
dat$Tribal_religious <- as.factor(dat$Tribal_religious)

# Benefit
dat <- dat |> 
  mutate(Benefit = case_when(Q16 == 1 ~ "Social order",
                             Q16 == 2 ~ "Biologial weak",
                             Q16 == 3 ~ "Poor"))
dat$Benefit <- as.factor(dat$Benefit)

# gender
dat <- dat |> 
  mutate(Gender = case_when(Gender == 1 ~ "Male",
                            Gender == 2 ~ "Female"))
dat$Gender <- as.factor(dat$Gender)

# urban/rural
dat <- dat |> 
  mutate(Urban_rural = case_when(D8 == 1 ~ "Urban",
                                 D8 == 2 ~ "Rural"))
dat$Urban_rural <- as.factor(dat$Urban_rural)

# sect
dat <- dat |> 
  mutate(Sect = case_when(D6 == 1 ~ "Sunni",
                          D6 == 2 ~ "Shia",
                          D6 == 3 ~ "Christian",
                          D6 == 4 ~ "Kurd",
                          D6 == 5 ~ "others"))
dat$Sect <- as.factor(dat$Sect)

# patronage_religious_daily_necessities
dat <- dat |> 
  mutate(ptr_rel_daily = case_when(Q4 >= 4 ~ "High",
                                   Q4 <= 2 ~ "Low",
                                   Q4 == 3 ~ "other"))
dat$ptr_rel_daily <- as.factor(dat$ptr_rel_daily)

# patronage_religious_job
dat <- dat |> 
  mutate(ptr_rel_job = case_when(Q5 >= 4 ~ "High",
                                 Q5 <= 2 ~ "Low",
                                 Q5 == 3 ~ "other"))
dat$ptr_rel_job <- as.factor(dat$ptr_rel_job)

# patronage_Tribe_daily_necessities
dat <- dat |> 
  mutate(ptr_trb_daily = case_when(Q6 >= 4 ~ "High",
                                   Q6 <= 2 ~ "Low",
                                   Q6 == 3 ~ "other"))
dat$ptr_trb_daily <- as.factor(dat$ptr_trb_daily)

# patronage_Tribe_job
dat <- dat |> 
  mutate(ptr_trb_job = case_when(Q7 >= 4 ~ "High",
                                 Q7 <= 2 ~ "Low",
                                 Q7 == 3 ~ "other"))
dat$ptr_trb_job <- as.factor(dat$ptr_trb_job)


# Liberal_high_education
dat <- dat |> 
  mutate(Education = case_when(D3 <= 2 ~ "Low",
                               D3 == 3 | D3 == 4 ~ "Middle",
                               D3 >= 5 ~ "High"))
dat$Education <- as.factor(dat$Education)

dat <- dat |>
  mutate(Liberal_high_education = ifelse(Liberalism == "High" & Education == "High", "High-educated_liberal", "Other"))
dat$Liberal_high_education <- as.factor(dat$Liberal_high_education)

# Shia parties supporter
dat <- dat |> 
  mutate(Shia_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Shia_party_supporter = ifelse(Shia_party_supporter >= 7, "Shia Party Supporter", "others"))
dat$Shia_party_supporter <- as.factor(dat$Shia_party_supporter)

# trust tribal chief
dat <- dat |> 
  mutate(trust_tribalchief = case_when(Q6_9 <= 2 ~ "High",
                                       Q6_9 == 4 | Q6_9 == 5 ~ "Low",
                                       Q6_9 == 3 | Q6_9 == 6 ~ "other"))
dat$trust_tribalchief <- as.factor(dat$trust_tribalchief)

# trust religious leaders
dat <- dat |> 
  mutate(trust_religiousleader = case_when(Q6_8 <= 2 ~ "High",
                                           Q6_8 == 4 | Q6_8 == 5 ~ "Low",
                                           Q6_8 == 3 | Q6_8 == 6 ~ "other"))
dat$trust_religiousleader <- as.factor(dat$trust_religiousleader)

# region south
dat <- dat |> 
  mutate(South = ifelse(D7 >= 12, "South", "other"))
dat$South <- as.factor(dat$South)

# rural south
dat <- dat |>
  mutate(rural_south = ifelse(South == "South" & Urban_rural == "Rural", "Rural_south", "Other"))
dat$rural_south <- as.factor(dat$rural_south)

# region Arabat
dat <- dat |> 
  mutate(Atabat = ifelse(D7 == 12 | D7 == 13, "Atabat", "other"))
dat$Atabat <- as.factor(dat$Atabat)

# region west
dat <- dat |> 
  mutate(West = ifelse(D7 >= 4 & D7 <= 8, "West", "other"))
dat$West <- as.factor(dat$West)

# region KRG
dat <- dat |> 
  mutate(KRG = ifelse(D7 <= 3, "KRG", "other"))
dat$KRG <- as.factor(dat$KRG)

# South Islamism
dat <- dat |>
  mutate(South_Islamism = ifelse(Islamism == "High" & South == "South", "South_Islamism", "Other"))
dat$South_Islamism <- as.factor(dat$South_Islamism)

# West Tribalism
dat <- dat |>
  mutate(West_tribalism = ifelse(Tribalism == "High" & West == "West", "West_Tribalism", "Other"))
dat$West_tribalism <- as.factor(dat$West_tribalism)

# South Tribalism
dat <- dat |>
  mutate(South_tribalism = ifelse(Tribalism == "High" & South == "South", "South_Tribalism", "Other"))
dat$South_tribalism <- as.factor(dat$South_tribalism)

# rural Islamism
dat <- dat |>
  mutate(rural_Islamism = ifelse(Islamism == "High" & Urban_rural == "Rural", "Rural_Islamism", "Other"))
dat$rural_Islamism <- as.factor(dat$rural_Islamism)

# rural Tribalism
dat <- dat |>
  mutate(rural_tribalism = ifelse(Tribalism == "High" & Urban_rural == "Rural", "Rural_tribalism", "Other"))
dat$rural_tribalism <- as.factor(dat$rural_tribalism)

# task
dat$task <- as.factor(dat$task)

```

# conjoint setting
```{r}
f1 <- selected ~ General_description + Religious_customary + Support_criticism + Internal_external
str(dat)

dat$General_description <- factor(dat$General_description, 
                          levels = c("De fact legalization","Right to choose"))
dat$Religious_customary <- factor(dat$Religious_customary, 
                          levels = c("Religious laws", "Tribal customs"))
dat$Support_criticism <- factor(dat$Support_criticism, 
                          levels = c("Legal fragmentaion", "Protect poor women", "Social order", "Violation to women's right"))
dat$Internal_external <- factor(dat$Internal_external, 
                          levels = c("Domestic criticism", "International criticism", "Evaluation"))
```

# whole model
```{r, fig.height = 5, dpi = 300}
# MM
whole_mm <- cj(dat, f1, 
               id = id,
               estimate = "mm", 
               h0 = 0.5)
whole_mm
plot(whole_mm, vline = 0.5, 
     size = 2)
write_html(whole_mm, file = "result/whole_mm.html")

# amce
whole_amce <- cj(data = dat, f1,
                 id = ~ respondentIndex,
                 estimate = "amce")
whole_amce
whole_model_amce <- plot(whole_amce, vline = 0, size = 2)
whole_model_amce

ggsave(filename = "result/whole_amce.png",
       plot = whole_model_amce,
       width = 8.5,
       height = 4.5,
       units = "in", 
       dpi = 600)

write_html(whole_amce, file = "result/whole_amce.html")

# plot estimation difference
whole_mm$Estimate <- "MM"
whole_amce$Estimate <- "AMCE"
mm_amce <- plot(rbind(whole_mm, whole_amce), feature_headers = FALSE, size = 2)+
  ggplot2::facet_wrap(~ Estimate, ncol = 2L)
mm_amce

ggsave(filename = "result/whole_mm_amce.png",
       plot = mm_amce,
       width = 8.5,
       height = 4.5,
       units = "in", 
       dpi = 600)
```

# Gender
```{r, fig.height = 5, dpi = 300}
# amce
Gender_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Gender)
Gender_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Gender)
plot(rbind(Gender_amce, Gender_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$Gender <- relevel(dat$Gender, "Male")
Gender_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ Gender)
Gender_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Gender)
plot(rbind(Gender_mm, Gender_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# plot gender
```{r, fig.height = 5, dpi = 300}
plot(Gender_amce, group = "Gender")

gender_mm <- plot(Gender_mm, group = "Gender", size = 2)
gender_mm
ggsave(filename = "result/gender_mm.png",
       plot = gender_mm,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

plot(Gender_amce, size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(Gender_mm, size = 2, vline = 0.5) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Gender_amce_diff, size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(Gender_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 2L)
```

# OLS result of Gender
```{r, fig.height = 5, dpi = 300}

write_html(Gender_amce, file = "result/Gender_ame.html")
write_html(Gender_mm, file = "result/Gender_mm.html")
write_html(Gender_amce_diff, file = "result/Gender_ame_diff.html")
write_html(Gender_mm_diff, file = "result/Gender_mm_diff.html")

# anova
cj_anova(data = dat, 
         selected ~ General_description + Religious_customary + Support_criticism + Internal_external, by = ~ Gender)
```

# Childe_marriage
```{r, fig.height = 5, dpi = 300}
table(dat$Childe_marriage)
# amce
Childe_marriage_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Childe_marriage)
Childe_marriage_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Childe_marriage)
plot(rbind(Childe_marriage_amce, Childe_marriage_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$Childe_marriage <- relevel(dat$Childe_marriage, "Disagree")
Childe_marriage_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ Childe_marriage)
Childe_marriage_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Childe_marriage)
plot(rbind(Childe_marriage_mm, Childe_marriage_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

plot(Childe_marriage_amce, group = "Childe_marriage")

child_marriage_mm <- plot(Childe_marriage_mm, group = "Childe_marriage", size = 2)
child_marriage_mm
ggsave(filename = "result/child_marriage_mm.png",
       plot = child_marriage_mm,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

```

# OLS result of childe marriage
```{r, fig.height = 5, dpi = 300}
write_html(Childe_marriage_amce, file = "result/Childe_marriage_ame.html")
write_html(Childe_marriage_mm, file = "result/Childe_marriage_mm.html")
write_html(Childe_marriage_amce_diff, file = "result/Childe_marriage_ame_diff.html")
write_html(Childe_marriage_mm_diff, file = "result/Childe_marriager_mm_diff.html")

# anova
cj_anova(data = dat, 
         selected ~ General_description + Religious_customary + Support_criticism + Internal_external, by = ~ Childe_marriage)
```

# H1a: Islamism
```{r, fig.height = 5, dpi = 300}
# amce
Islam_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Islamism)
Islam_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Islamism)
plot(rbind(Islam_amce, Islam_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$Islamism <- relevel(dat$Islamism, "Low")
Islam_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ Islamism)
Islam_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Islamism)
plot(rbind(Islam_mm, Islam_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

Islamism <- plot(Islam_mm, group = "Islamism", size = 2, vline = 0)
Islamism
ggsave(filename = "result/Islamism_mm.png",
       plot = Islamism,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

Islamism_diff <- plot(Islam_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 2L)
Islamism_diff
ggsave(filename = "result/Islamism_diff.png",
       plot = Islamism_diff,
       width = 9,
       height = 4,
       units = "in", 
       dpi = 600)
```

# OLS result of Islamism
```{r, fig.height = 5, dpi = 300}
write_html(Islam_amce, file = "result/Islam_ame.html")
write_html(Islam_mm, file = "result/Islame_mm.html")
write_html(Islam_amce_diff, file = "result/Islam_ame_diff.html")
write_html(Islam_mm_diff, file = "result/Islam_mm_diff.html")

# anova
cj_anova(data = dat, 
         selected ~ General_description + Religious_customary + Support_criticism + Internal_external, by = ~ Islamism)
```

# H1b: Shia parties supporter
```{r, fig.height = 5, dpi = 300}
# amce
Shia_party_supporter_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ Shia_party_supporter)
Shia_party_supporter_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ Shia_party_supporter)
plot(rbind(Shia_party_supporter_amce, Shia_party_supporter_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
Shia_party_supporter_mm <- cj(dat,
                       selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                       id = id,
                       estimate = "mm",
                       by = ~ Shia_party_supporter)
Shia_party_supporter_mm_diff <- cj(data = dat,
                            selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                            id = id,
                            estimate = "mm_diff",
                            by = ~ Shia_party_supporter)
plot(rbind(Shia_party_supporter_mm, Shia_party_supporter_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

party <- plot(Shia_party_supporter_mm, group = "Shia_party_supporter", size = 2)
party
ggsave(filename = "result/party_mm.png",
       plot = party,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

party_diff <- plot(Shia_party_supporter_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 2L)
party_diff
ggsave(filename = "result/party_diff.png",
       plot = party_diff,
       width = 9,
       height = 4,
       units = "in", 
       dpi = 600)
```

# OLS result of Shia_party_supporter
```{r, fig.height = 5, dpi = 300}
write_html(Shia_party_supporter_amce, file = "result/Shia_party_supporter_ame.html")
write_html(Shia_party_supporter_mm, file = "result/Shia_party_supporter_mm.html")
write_html(Shia_party_supporter_amce_diff, file = "result/Shia_party_supporter_ame_diff.html")
write_html(Shia_party_supporter_mm_diff, file = "result/Shia_party_supporter_mm_diff.html")

# anova
cj_anova(data = dat, 
         selected ~ General_description + Religious_customary + Support_criticism + Internal_external, by = ~ Shia_party_supporter)
```

# H1c: Patronage: Religioun and daily necessities
```{r, fig.height = 5, dpi = 300}
# amce
ptr_rel_daily_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ ptr_rel_daily)
ptr_rel_daily_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ ptr_rel_daily)
plot(rbind(ptr_rel_daily_amce, ptr_rel_daily_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$ptr_rel_daily <- relevel(dat$ptr_rel_daily, "Low")
ptr_rel_daily_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ ptr_rel_daily)
ptr_rel_daily_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ ptr_rel_daily)
plot(rbind(ptr_rel_daily_mm, ptr_rel_daily_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

ptr_rel <- plot(ptr_rel_daily_mm, group = "ptr_rel_daily", size = 2)
ptr_rel
ggsave(filename = "result/Patronage_rel_mm.png",
       plot = ptr_rel,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

ptr_rel_diff <- plot(ptr_rel_daily_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 2L)
ptr_rel_diff
ggsave(filename = "result/Patronage_rel_diff.png",
       plot = ptr_rel_diff,
       width = 9,
       height = 4,
       units = "in", 
       dpi = 600)
```

# OLS result of Patronage: Religioun and daily necessities
```{r, fig.height = 5, dpi = 300}
write_html(ptr_rel_daily_amce, file = "result/ptr_rel_daily_ame.html")
write_html(ptr_rel_daily_mm, file = "result/ptr_rel_daily_mm.html")
write_html(ptr_rel_daily_amce_diff, file = "result/ptr_rel_daily_ame_diff.html")
write_html(ptr_rel_daily_mm_diff, file = "result/ptr_rel_daily_mm_diff.html")

# anova
cj_anova(data = dat, 
         selected ~ General_description + Religious_customary + Support_criticism + Internal_external, by = ~ ptr_rel_daily)
```

# H2a: Tribalism
```{r, fig.height = 5, dpi = 300}
# amce
Tribalism_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Tribalism)
Tribalism_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Tribalism)
plot(rbind(Tribalism_amce, Tribalism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$Tribalism <- relevel(dat$Tribalism, "Low")
Tribalism_mm <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "mm",
                 by = ~ Tribalism)
Tribalism_mm_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "mm_diff",
                      by = ~ Tribalism)
plot(rbind(Tribalism_mm, Tribalism_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

tribal <- plot(Tribalism_mm, group = "Tribalism", size = 2)
tribal
ggsave(filename = "result/Tribalism_mm.png",
       plot = tribal,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

tribal_diff <- plot(Tribalism_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 2L)
tribal_diff
ggsave(filename = "result/Tribalism_diff.png",
       plot = tribal_diff,
       width = 9,
       height = 4,
       units = "in", 
       dpi = 600)
```

# OLS result of Tribalism
```{r, fig.height = 5, dpi = 300}
write_html(Tribalism_amce, file = "result/Tribalism_ame.html")
write_html(Tribalism_mm, file = "result/Tribalism_mm.html")
write_html(Tribalism_amce_diff, file = "result/Tribalism_diff.html")
write_html(Tribalism_mm_diff, file = "result/Tribalism_mm_diff.html")

# anova
cj_anova(data = dat, 
         selected ~ General_description + Religious_customary + Support_criticism + Internal_external, by = ~ Tribalism)
```

# H2b: Patronage: Triba and daily necessities
```{r, fig.height = 5, dpi = 300}
table(dat$ptr_trb_daily)
# amce
ptr_trb_daily_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ ptr_trb_daily)
ptr_trb_daily_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ ptr_trb_daily)
plot(rbind(ptr_trb_daily_amce, ptr_trb_daily_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$ptr_trb_daily <- relevel(dat$ptr_trb_daily, "Low")
ptr_trb_daily_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ ptr_trb_daily)
ptr_trb_daily_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ ptr_trb_daily)
plot(rbind(ptr_trb_daily_mm, ptr_trb_daily_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

ptr_trb <- plot(ptr_trb_daily_mm, group = "ptr_trb_daily", size = 2)
ptr_trb
ggsave(filename = "result/Patronage_trb_mm.png",
       plot = ptr_trb,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)

ptr_trb_diff <- plot(ptr_trb_daily_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~ BY, ncol = 2L)
ptr_trb_diff
ggsave(filename = "result/Patronage_trb_diff.png",
       plot = ptr_trb_diff,
       width = 9,
       height = 4,
       units = "in", 
       dpi = 600)
```

# OLS result of Patronage: Triba and daily necessities
```{r, fig.height = 5, dpi = 300}
write_html(ptr_trb_daily_amce, file = "result/ptr_trb_daily_ame.html")
write_html(ptr_trb_daily_mm, file = "result/ptr_trb_daily_mm.html")
write_html(ptr_trb_daily_amce_diff, file = "result/ptr_trb_daily_diff.html")
write_html(ptr_trb_daily_mm_diff, file = "result/ptr_trb_daily_mm_diff.html")

# anova
cj_anova(data = dat, 
         selected ~ General_description + Religious_customary + Support_criticism + Internal_external, by = ~ ptr_trb_daily)
```

# frequency and propotion
```{r}
plot(cj_freqs(dat,
              f1,
              id = ~ respondentIndex))
```






# Rural Tribalism
```{r, fig.height = 5, dpi = 300}
# amce
rural_tribalism_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ rural_tribalism)
rural_tribalism_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ rural_tribalism)
plot(rbind(rural_tribalism_amce, rural_tribalism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
rural_tribalism_mm <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "mm",
                        by = ~ rural_tribalism)
rural_tribalism_mm_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "mm_diff",
                              by = ~ rural_tribalism)
plot(rbind(rural_tribalism_mm, rural_tribalism_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# task
```{r, fig.height = 7, dpi = 300}
# amce
task_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ task)
task_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ task)
plot(task_amce, size = 2, vline = 0) +
  ggplot2::facet_wrap(~ BY, ncol = 3L)
plot(rbind(task_amce, task_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Patronage: Triba and job
```{r, fig.height = 7, dpi = 300}
table(dat$ptr_trb_job)
# amce
ptr_trb_job_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ ptr_trb_job)
ptr_trb_job_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ ptr_trb_job)
plot(rbind(ptr_trb_daily_amce, ptr_trb_daily_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$ptr_trb_job <- relevel(dat$ptr_trb_job, "Low")
ptr_trb_job_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ ptr_trb_job)
ptr_trb_job_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ ptr_trb_job)
plot(rbind(ptr_trb_job_mm, ptr_trb_job_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# region West
```{r, fig.height = 5, dpi = 300}
# amce
West_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ West)
West_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ West)
plot(rbind(West_amce, West_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
West_mm <- cj(dat,
                       selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                       id = id,
                       estimate = "mm",
                       by = ~ West)
West_mm_diff <- cj(data = dat,
                            selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                            id = id,
                            estimate = "mm_diff",
                            by = ~ West)
plot(rbind(West_mm, West_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# region KRG
```{r, fig.height = 5, dpi = 300}
# amce
KRG_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ KRG)
KRG_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ KRG)
plot(rbind(KRG_amce, KRG_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
KRG_mm <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "mm",
                        by = ~ KRG)
KRG_mm_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "mm_diff",
                              by = ~ KRG)
plot(rbind(KRG_mm, KRG_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# South Tribalism
```{r, fig.height = 5, dpi = 300}
# amce
South_tribalism_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ South_tribalism)
South_tribalism_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ South_tribalism)
plot(rbind(South_tribalism_amce, South_tribalism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# MM
South_tribalism_mm <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "mm",
                        by = ~ South_tribalism)
South_tribalism_mm_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "mm_diff",
                              by = ~ South_tribalism)
plot(rbind(South_tribalism_mm, South_tribalism_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# Liberalism
```{r, fig.height = 7, dpi = 300}
# amce
Liberalism_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Liberalism)
Liberalism_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Liberalism)
plot(rbind(Liberalism_amce, Liberalism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$Liberalism <- relevel(dat$Liberalism, "Low")
Liberalism_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ Liberalism)
Liberalism_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Liberalism)
plot(rbind(Liberalism_mm, Liberalism_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```


# Tribal_religious
```{r, fig.height = 5, dpi = 300}
# amce
Tribal_religious_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Tribal_religious)
Tribal_religious_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Tribal_religious)
plot(rbind(Tribal_religious_amce, Tribal_religious_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Tribal_religious_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ Tribal_religious)
Tribal_religious_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Tribal_religious)
plot(rbind(Tribal_religious_mm, Tribal_religious_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Benefit
```{r, fig.height = 7, dpi = 300}
# amce
Benefit_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Benefit)
Benefit_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Benefit)
plot(rbind(Benefit_amce, Benefit_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$Benefit <- relevel(dat$Benefit, "Social order")
Benefit_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ Benefit)
Benefit_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Benefit)
plot(rbind(Benefit_mm, Benefit_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Urban_rural
```{r, fig.height = 5, dpi = 300}
# amce
Urban_rural_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Urban_rural)
Urban_rural_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Urban_rural)
plot(rbind(Urban_rural_amce, Urban_rural_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
Urban_rural_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ Urban_rural)
Urban_rural_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Urban_rural)
plot(rbind(Urban_rural_mm, Urban_rural_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```

# Sect
```{r, fig.height = 7, dpi = 300}
# amce
Sect_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Sect)
Sect_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Sect)
plot(rbind(Sect_amce, Sect_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

```

# Patronage: Religioun and job
```{r, fig.height = 7, dpi = 300}
# amce
ptr_rel_job_amce <- cj(data = dat,
                 selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ ptr_rel_job)
ptr_rel_job_amce_diff <- cj(data = dat,
                      selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ ptr_rel_job)
plot(rbind(ptr_rel_job_amce, ptr_rel_job_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)

# mm
dat$ptr_rel_job <- relevel(dat$ptr_rel_job, "Low")
ptr_rel_job_mm <- cj(dat,
               selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
               id = id,
               estimate = "mm",
               by = ~ ptr_rel_job)
ptr_rel_job_mm_diff <- cj(data = dat,
                    selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ ptr_rel_job)
plot(rbind(ptr_rel_job_mm, ptr_rel_job_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~ BY, ncol = 3L)
```


# Liberalism and high education
```{r, fig.height = 5, dpi = 300}
# amce
liberal_hiedu_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ Liberal_high_education)
liberal_hiedu_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ Liberal_high_education)
plot(rbind(liberal_hiedu_amce, liberal_hiedu_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
dat$Liberal_high_education <- relevel(dat$Liberal_high_education, "Other")

liberal_hiedu_mm <- cj(dat,
                       selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                       id = id,
                       estimate = "mm",
                       by = ~ Liberal_high_education)
liberal_hiedu_mm_diff <- cj(data = dat,
                            selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                            id = id,
                            estimate = "mm_diff",
                            by = ~ Liberal_high_education)
plot(rbind(liberal_hiedu_mm, liberal_hiedu_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```


# Trust tribal chief
```{r, fig.height = 7, dpi = 300}
# amce
trust_tribalchief_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ trust_tribalchief)
trust_tribalchief_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ trust_tribalchief)
plot(rbind(trust_tribalchief_amce, trust_tribalchief_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
trust_tribalchief_mm <- cj(dat,
                       selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                       id = id,
                       estimate = "mm",
                       by = ~ trust_tribalchief)
trust_tribalchief_mm_diff <- cj(data = dat,
                            selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                            id = id,
                            estimate = "mm_diff",
                            by = ~ trust_tribalchief)
plot(rbind(trust_tribalchief_mm, trust_tribalchief_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# Trust tribal religious leaders
```{r, fig.height = 7, dpi = 300}
# amce
trust_religiousleader_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ trust_religiousleader)
trust_religiousleader_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ trust_religiousleader)
plot(rbind(trust_religiousleader_amce, trust_religiousleader_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

```

# Region south
```{r, fig.height = 5, dpi = 300}
# amce
South_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ South)
South_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ South)
plot(rbind(South_amce, South_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
Shia_party_supporter_mm <- cj(dat,
                       selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                       id = id,
                       estimate = "mm",
                       by = ~ Shia_party_supporter)
Shia_party_supporter_mm_diff <- cj(data = dat,
                            selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                            id = id,
                            estimate = "mm_diff",
                            by = ~ Shia_party_supporter)
plot(rbind(Shia_party_supporter_mm, Shia_party_supporter_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# region rural south
```{r, fig.height = 5, dpi = 300}
# amce
rural_south_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ rural_south)
rural_south_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ rural_south)
plot(rbind(rural_south_amce, rural_south_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# region Atabat
```{r, fig.height = 5, dpi = 300}
# amce
Atabat_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ Atabat)
Atabat_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ Atabat)
plot(rbind(Atabat_amce, Atabat_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```


# South Islamism
```{r, fig.height = 5, dpi = 300}
# amce
South_Islamism_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ South_Islamism)
South_Islamism_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ South_Islamism)
plot(rbind(South_Islamism_amce, South_Islamism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```


# West Tribalism
```{r, fig.height = 5, dpi = 300}
# amce
West_tribalism_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ West_tribalism)
West_tribalism_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ West_tribalism)
plot(rbind(West_tribalism_amce, West_tribalism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# Rural Islamism
```{r, fig.height = 5, dpi = 300}
# amce
rural_Islamism_amce <- cj(data = dat,
                        selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ rural_Islamism)
rural_Islamism_amce_diff <- cj(data = dat,
                              selected ~ General_description + Religious_customary + Support_criticism + Internal_external,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ rural_Islamism)
plot(rbind(rural_Islamism_amce, rural_Islamism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```


